scRNAseq analysis for Deng et al. (2014)
Jose Alejandro Romero Herrera
2021-12-17
1 Introduction
scRNAseq data analysis of Deng et al. (2014). This study follows the early embryonic development of mouse cells from zygote to late blastocists (and some adult cells).
Mouse Early Development https://doi.org/10.1093/molehr/gaw002
2 Load libraries
library(rhdf5)
library(tidyverse)
library(Seurat)
library(SeuratWrappers)
library(patchwork)
library(biomaRt)
library(slingshot)
library(gprofiler2)3 Load dataset
Data is a matrix of cells in columns and genes in rows. We specify that the row names (our genes) are in the first column.
raw_data <- read.delim("../Data/example_data/GSE45719_cts.txt", sep = "\t", row.names = 1)
head(raw_data, n = 5)Here is the metadata. The column names of raw_data are the same as the row names of metadata.
metadata <- read.delim("../Data/example_data/GSE45719_metadata.tsv", header = T, sep = "\t", quote = "", stringsAsFactors = F, row.names = 1)
head(metadata, n = 5)For the purposes of this workshop, the metadata of this project has been simplified. The most interesting variable is Cell_type, which shows the developmental stage of our cells. Cell type is not ordered by its differentiation stage, but we can do that now:
metadata$Cell_type <- factor(metadata$Cell_type,
levels = c("MII Oocyte","Zygote",
"Early 2-cell","Mid 2-cell","Late 2-cell",
"4-cell","8-cell","16-cell",
"Early blastocyst","Mid blastocyst","Late blastocyst",
"Fibroblast","Adult"))There are some cell types that we are not really interested in, and will be removed in the next step below.
3.1 Seurat object
Now that we have our data loaded, we can create a Seurat Object, from which we will perform our analysis. The function CreateSeuratObject() takes as arguments our count matrix, metadata and several options for filtering the data:
min.cellswill keep features (genes) that are detected in at least this many cells.min.featureswill keep cells with at least this many features detected.
This will prevent from keeping cells and genes with an immense majority of 0’s as values. In addition, we are removing some cell types we do not want to analyse.
raw_ann <- CreateSeuratObject(counts = raw_data, meta.data = metadata, min.cells = 3, min.features = 200,)
raw_ann <- subset(raw_ann, cells = colnames(raw_ann)[!raw_ann$Cell_type %in% c("MII Oocyte","Fibroblast","Adult")])4 Quality Control
The Seurat object initialization step above only considered cells that expressed at least 300 genes and genes detected in at least 3 cells. Here is how many cells and genes we start with:
print(paste0("Before filtering: ", dim(raw_ann)[2], " cells ", dim(raw_ann)[1], " genes"))## [1] "Before filtering: 286 cells 20566 genes"Additionally, we would like to exclude cells that are damaged. A common metric to judge this (although by no means the only one) is the relative expression of mitochondrially derived genes. When the cells apoptose due to stress, their mitochondria becomes leaky and there is widespread RNA degradation. Thus a relative enrichment of mitochondrially derived genes can be a tell-tale sign of cell stress. Here, we compute the proportion of transcripts that are of mitochondrial origin for every cell (percent.mito), and visualize its distribution as a violin plot. We also use the GenePlot() function to observe how percent.mito correlates with other metrics.
raw_ann[['percent.mito']] <- PercentageFeatureSet(raw_ann, pattern = "^mt-")Now, sometimes, mitochondrial genes are a bit tricky to find, specially if your genes are not gene names, but gene IDs. You might have already a collection of gene names you want to use, but it is not always the case. Therefore, it might very useful to have some extra annotation on our genes which will help to select mitochondrial genes (or ERCC genes and ribosomal genes). Since this dataset was aligned using the mouse genome version mm9, we will use the mm9 annotation from Biomart. The package biomaRt will do the job!
ensembl67=useMart(host='http://feb2014.archive.ensembl.org/', # select latest version of mm9 genome annotation
biomart='ENSEMBL_MART_ENSEMBL', dataset = "mmusculus_gene_ensembl")
# we get a list of annotations we would like to fetch
mm9.gene.annotations <- biomaRt::getBM(mart = ensembl67, attributes=c("ensembl_gene_id", "external_gene_id", "description", "chromosome_name"))
head(mm9.gene.annotations)We select that mitochondrial genes are those with the chromosome_name annotation equal to “MT”.
mt_genes <- mm9.gene.annotations %>% filter(chromosome_name == "MT") %>% pull(external_gene_id)
mt_genes <- mt_genes %in% rownames(raw_ann)
raw_ann[['percent.mito']] <- PercentageFeatureSet(raw_ann, features = mt_genes)If your data contains spike-ins (ERCC genes), you can also compute the proportion of transcripts belonging to them. Furthermore, it might be useful to calculate the percentage of reads aligned to ribosomal genes, since it has been shown that they can skew the data due to their high variability.
raw_ann[['percent.ercc']] <- PercentageFeatureSet(raw_ann, pattern = "^ERCC-")
raw_ann[['percent.ribo']] <- PercentageFeatureSet(raw_ann, pattern = "^Rp[ls]")Nonetheless, in this experiment there aren’t ERCC genes or mitochondrial genes:
sum(raw_ann$percent.ercc)## [1] 0sum(raw_ann$percent.mito)## [1] 04.1 Visualizations
It is extremely useful to visualize QC measurements calculated so far. Violin plots (fancy boxplots) are a great way to check the distribution of values of all our QC measurements. When initializing the Seurat Object, Seurat calculates also the number of genes detected and the total library size per cell.
VlnPlot(raw_ann,
features = c("nFeature_RNA", "nCount_RNA", "percent.ribo"),
ncol = 4, group.by = "Cell_type")
It seems that the late 2-cells have a very high total number of reads. We should not filter them out, there are very few! In addition, our percentage of ribosomal counts is also quite low (maximum is ~5%).
We can check the relationship between library size and number of genes detected:
FeatureScatter(raw_ann, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "Cell_type")
In other experiment with more cells, we might consider removing cells with an unexpectedly high number of detected genes and reads. These cells might be doublets, that is, two cells that were sequenced together!
4.2 Filtering
Some automatic filtering can be made using quantiles! Using the extreme upper and lower quantiles (0.99 and 0.01) we can make sure that outliers are removed. In this case, we will remove cells with very low library size and genes detected (lower than quantile 0.01).
feature_min <- quantile(raw_ann$nFeature_RNA, probs = 0.01)
count_min <- quantile(raw_ann$nCount_RNA, probs = 0.01)We can subset our dataset using the function subset().
adata <- subset(raw_ann, subset =
nFeature_RNA > feature_min &
nCount_RNA > count_min)
rm(raw_ann) # we remove the initial unfiltered dataset to reduce computational resources, this is not necessary!Finally, this is how many cells and genes we have after filtering:
print(paste0("After filtering: ", dim(adata)[2], " cells ", dim(adata)[1], " genes"))## [1] "After filtering: 280 cells 20566 genes"And this is how our filtered data looks like:
VlnPlot(adata,
features = c("nFeature_RNA", "nCount_RNA", "percent.ribo"),
ncol = 4, group.by = "Cell_type")
5 Normalization
Now that the data is filtered, we can proceed to normalize our count matrix. Seurat normalizes the gene expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. There have been many methods to normalize the data, but this is the simplest and the most intuitive. The division by total expression is done to change all expression counts to a relative measure, since experience has suggested that technical factors (e.g. capture rate, efficiency of reverse transcription) are largely responsible for the variation in the number of molecules per cell, although genuine biological factors (e.g. cell cycle stage, cell size) also play a smaller, but non-negligible role. The log-transformation is a commonly used transformation that has many desirable properties, such as variance stabilization (can you think of others?).
adata <- NormalizeData(adata)
adata <- FindVariableFeatures(adata, selection.method = "vst", nfeatures = 2000)6 Identification of Variable Genes
Identify most variable genes and label top 5 most highly variable.
top10 <- head(VariableFeatures(adata), 5)
LabelPoints(plot = VariableFeaturePlot(adata), points = top10, repel = T)
7 Scale
adata <- ScaleData(adata)8 PCA
adata <- RunPCA(adata)pcalod_1 <- VizDimLoadings(object = adata, dims = 1) + theme(axis.text.y = element_text(size = 8))
pcalod_2 <- VizDimLoadings(object = adata, dims = 2) + theme(axis.text.y = element_text(size = 8))CombinePlots(plots = list(pcalod_1, pcalod_2), ncol = 2)
Visualization only show negative loadings cause there are many more than positives. We can see a balanced plot using balance = TRUE
pcalod_1 <- VizDimLoadings(object = adata, dims = 1, balanced = TRUE) + theme(axis.text.y = element_text(size = 8))
pcalod_2 <- VizDimLoadings(object = adata, dims = 2, balanced = TRUE) + theme(axis.text.y = element_text(size = 8))CombinePlots(plots = list(pcalod_1, pcalod_2), ncol = 2)
9 Clustering
Jackstraw method and plot to determine proper number of dimensions. Including elbow plot
adata <- JackStraw(adata, num.replicate = 100)
adata <- ScoreJackStraw(adata, dims = 1:20)
JackStrawPlot(adata, dims = 1:20)
ElbowPlot(adata)
We could include 6 PCs, but it does not hurt to use more, there are very few cells in this experiment.
adata <- FindNeighbors(adata, dims = 1:20)
adata <- FindClusters(adata)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 280
## Number of edges: 6771
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7890
## Number of communities: 7
## Elapsed time: 0 seconds10 Visualization
adata <- RunTSNE(adata)
adata <- RunUMAP(adata, dims = 1:20)10.1 PCA
DimPlot(adata, reduction = "pca", group.by = c('Cell_type')) +
DimPlot(adata, reduction = "pca", group.by = c('seurat_clusters'))
10.2 UMAP
p1 <- DimPlot(adata, reduction = "umap", group.by = 'Cell_type')
p2 <- DimPlot(adata, reduction = "umap")
p1 + p2
10.3 TSNE
p3 <- DimPlot(adata, reduction = "tsne", group.by = 'Cell_type')
p4 <- DimPlot(adata, reduction = "tsne")
p3 + p4
11 Differential Expression Analysis
Differential testing works a bit different than in bulk RNA-Seq analysis. Usually, your experiment will include dozens or hundreds of cells per cluster/condition. We can make use of Wilcox tests and multiple testing correction to identify statistically significant genes using the FindAllMarkers() function. This function will gather all “markers” for each of your Ident variable values (your conditions or identified clusters).
markers <- FindAllMarkers(adata,logfc.threshold = 0.5, only.pos = T)On the other hand, if you want to make a specific comparison, you can use the FindMarkers() function:
cluster2.markers <- FindMarkers(adata, ident.1 = 2, min.pct = 0.5, only.pos = T) # Only cluster 2 markersEven compare specific clusters against others by selecting ident.1 and ident.2:
cluster5.markers <- FindMarkers(adata, ident.1 = 5, ident.2 = c(0, 3), # Differences between cluster 5 and clusters 0 and 3
min.pct = 0.5, only.pos = T) 11.1 Visualizations
We can visualize all markers in a heatmap just like this! To go easy on the visualization, we only use the top 10 markers.
markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(adata, features = top10$gene) + NoLegend()
We can also plot the expression on our PCA, tSNE or UMAP. We only need to pass a vector of genes to the features argument of FeaturePlot()
FeaturePlot(adata, features = top10$gene[1:4], reduction = "pca")
Or as a violin plot, a ridge plot or a dot plot
VlnPlot(adata, features = top10$gene[1]) + RidgePlot(adata, features = top10$gene[1])
DotPlot(adata, features = top10$gene[1:10])
12 Functional analysis with gprofiler2
gost() function allows us to do functional profiling of gene lists, such as our differentially expressed genes. The function performs statistical enrichment analysis to find over-representation of terms from Gene Ontology, biological pathways like KEGG and Reactome, human disease annotations, etc. This is done by using hypergeometric tests that are corrected for multiple testing.
12.1 Single query
A standard input of the gost() function is a (named) list of gene identifiers. The list can consist of mixed types of identifiers (proteins, transcripts, microarray IDs, etc), SNP IDs, chromosomal intervals or functional term IDs.
The result is a named list where result is a data.frame with the enrichment analysis results and meta containing a named list with all the metadata for the query.
gostres <- gost(query = rownames(cluster2.markers),
organism = "mmusculus", ordered_query = FALSE,
multi_query = FALSE, significant = FALSE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.05, correction_method = "g_SCS",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
head(gostres$result)The result data.frame contains the following columns:
names(gostres$meta)## [1] "query_metadata" "result_metadata" "genes_metadata" "timestamp"
## [5] "version"12.2 Multiple queries
The function gost() also allows to perform enrichment on multiple input gene lists. Multiple queries are automatically detected if the input query is a list of vectors with gene identifiers and the results are combined into identical data.frame as in case of single query.
multi_gostres1 <- gost(query = list("chromX" = c("X:1000:1000000", "rs17396340",
"GO:0005005", "ENSG00000156103", "NLRP1"),
"chromY" = c("Y:1:10000000", "rs17396340",
"GO:0005005", "ENSG00000156103", "NLRP1")),
multi_query = FALSE)
head(multi_gostres1$result, 3)The column “query” in the result dataframe will now contain the corresponding name for the query. If no name is specified, then the query name is defined as the order of query with the prefix “query_.” Another option for multiple gene lists is setting the parameter multiquery = TRUE. Then the results from all of the input queries are grouped according to term IDs for better comparison.
multi_gostres2 <- gost(query = list("chromX" = c("X:1000:1000000", "rs17396340",
"GO:0005005", "ENSG00000156103", "NLRP1"),
"chromY" = c("Y:1:10000000", "rs17396340",
"GO:0005005", "ENSG00000156103", "NLRP1")),
multi_query = TRUE)
head(multi_gostres2$result, 3)12.3 Visualization
The enrichment results are visualized with a Manhattan-like-plot using the function gostplot() and the previously found gost results gostres:
#gostplot(gostres, capped = TRUE, interactive = TRUE)The function publish_gostplot() takes the static plot object as an input and enables to highlight a selection of interesting terms from the results with numbers and table of results. These can be set with parameter highlight_terms listing the term IDs in a vector or as a data.frame with column “term_id” such as a subset of the result dataframe.
First we create the static plot
p <- gostplot(gostres, capped = FALSE, interactive = FALSE)
p
Then we make it in high quality. We can add highlighted terms if we want with the highlight_terms argument.
pp <- publish_gostplot(p, highlight_terms = c("CORUM:3047"),
width = NA, height = NA, filename = NULL )
The gost results can also be visualized with a table. The publish_gosttable() function will create a nice-looking table with the result statistics for the highlight_terms from the result data.frame. The highlight_terms can be a vector of term IDs or a subset of the results.
publish_gosttable(gostres, highlight_terms = gostres$result[c(1:10),],
use_colors = TRUE,
show_columns = c("source", "term_name", "term_size", "intersection_size"),
filename = NULL) 
13 Pseudotime
We will use the package Slignshot. Slingshot was designed to model developmental trajectories in single-cell RNA sequencing data and serve as a component in an analysis pipeline after dimensionality reduction and clustering.
adata_SCE <- as.SingleCellExperiment(adata)
sde <- slingshot(adata_SCE, reducedDim = 'PCA')The using the slingshot package it is possible to plot our components and the calculated trajectory. Unfortunately, it does not work with ggplot.
cols <- colorRampPalette(viridis::viridis(3))
cols <- cols(50)[cut(sde$slingPseudotime_1, breaks=50)]
plot(reducedDim(adata_SCE, "PCA"), pch=16, col = cols) # adata_SCE contains, PC1 and PC2, and colored by pseudotime
lines(slingshot::SlingshotDataSet(sde), lwd = 2) # here we plot the trajectory
It is possible to create a FeaturePlot() with both pseudotime and the trajectory by extracting the information from our slingshot object and add it to our Seurat object, but it is a bit inconvenient. First, pseudotime can be extracted with the function slingPseudotime(). The trajectory curves are stored inside the slingshot object sde; in this case, it is easier since there is only one trajectory.
adata$pseudotime <- slingPseudotime(sde)
trajectory <- data.frame(SlingshotDataSet(sde)@curves[["Lineage1"]][["s"]])
head(trajectory)Now we can use Seurat functions to plot the pseudotime and the calculated trajectory!
FeaturePlot(adata, features = "pseudotime", reduction = "pca") + scale_color_viridis_c() + geom_path(data = trajectory, aes(x = PC_1, y = PC_2))
DimPlot(adata, group.by = "Cell_type", reduction = "pca")
We can create a plot that sorts our cells by pseudotime and check that it corresponds to their cell type
adata_SCE$pseudotime <- as.numeric(slingPseudotime(sde))
ggplot(as.data.frame(colData(adata_SCE)), aes(x = pseudotime,
y = Cell_type,
colour = Cell_type)) +
ggbeeswarm::geom_quasirandom(groupOnX = FALSE) +
theme_classic() +
xlab("Slingshot pseudotime") + ylab("Timepoint") +
ggtitle("Cells ordered by Slingshot pseudotime")
14 Session info
Finally, we create a session_info() table that will allow anyone to check what versions of R and packages are we using for reproducibility purposes.
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.1.0 (2021-05-18)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Copenhagen
## date 2021-12-17
## pandoc 2.14.0.1 @ /usr/local/bin/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.1.0)
## AnnotationDbi 1.56.2 2021-11-09 [1] Bioconductor
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
## backports 1.4.1 2021-12-13 [1] CRAN (R 4.1.0)
## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.1.0)
## Biobase * 2.54.0 2021-10-26 [1] Bioconductor
## BiocFileCache 2.2.0 2021-10-26 [1] Bioconductor
## BiocGenerics * 0.40.0 2021-10-26 [1] Bioconductor
## BiocManager 1.30.16 2021-06-15 [1] CRAN (R 4.1.0)
## biomaRt * 2.50.1 2021-11-21 [1] Bioconductor
## Biostrings 2.62.0 2021-10-26 [1] Bioconductor
## bit 4.0.4 2020-08-04 [1] CRAN (R 4.1.0)
## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.1.0)
## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.0)
## blob 1.2.2 2021-07-23 [1] CRAN (R 4.1.0)
## broom 0.7.10 2021-10-31 [1] CRAN (R 4.1.0)
## bslib 0.3.1 2021-10-06 [1] CRAN (R 4.1.0)
## cachem 1.0.6 2021-08-19 [1] CRAN (R 4.1.0)
## callr 3.7.0 2021-04-20 [1] CRAN (R 4.1.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
## cli 3.1.0 2021-10-27 [1] CRAN (R 4.1.0)
## cluster 2.1.2 2021-04-17 [1] CRAN (R 4.1.0)
## codetools 0.2-18 2020-11-04 [1] CRAN (R 4.1.0)
## colorspace 2.0-2 2021-06-24 [1] CRAN (R 4.1.0)
## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.1.0)
## crayon 1.4.2 2021-10-29 [1] CRAN (R 4.1.0)
## curl 4.3.2 2021-06-23 [1] CRAN (R 4.1.0)
## data.table 1.14.2 2021-09-27 [1] CRAN (R 4.1.0)
## DBI 1.1.1 2021-01-15 [1] CRAN (R 4.1.0)
## dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.1.0)
## DelayedArray 0.20.0 2021-10-26 [1] Bioconductor
## deldir 1.0-6 2021-10-23 [1] CRAN (R 4.1.0)
## desc 1.4.0 2021-09-28 [1] CRAN (R 4.1.0)
## devtools 2.4.3 2021-11-30 [1] CRAN (R 4.1.0)
## digest 0.6.29 2021-12-01 [1] CRAN (R 4.1.0)
## dplyr * 1.0.7 2021-06-18 [1] CRAN (R 4.1.0)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.1.0)
## fansi 0.5.0 2021-05-25 [1] CRAN (R 4.1.0)
## farver 2.1.0 2021-02-28 [1] CRAN (R 4.1.0)
## fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.1.0)
## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.1.0)
## fitdistrplus 1.1-6 2021-09-28 [1] CRAN (R 4.1.0)
## forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.1.0)
## fs 1.5.2 2021-12-08 [1] CRAN (R 4.1.0)
## future 1.23.0 2021-10-31 [1] CRAN (R 4.1.0)
## future.apply 1.8.1 2021-08-10 [1] CRAN (R 4.1.0)
## generics 0.1.1 2021-10-25 [1] CRAN (R 4.1.0)
## GenomeInfoDb * 1.30.0 2021-10-26 [1] Bioconductor
## GenomeInfoDbData 1.2.7 2021-11-16 [1] Bioconductor
## GenomicRanges * 1.46.1 2021-11-18 [1] Bioconductor
## ggbeeswarm 0.6.0 2017-08-07 [1] CRAN (R 4.1.0)
## ggplot2 * 3.3.5 2021-06-25 [1] CRAN (R 4.1.0)
## ggrepel 0.9.1 2021-01-15 [1] CRAN (R 4.1.0)
## ggridges 0.5.3 2021-01-08 [1] CRAN (R 4.1.0)
## globals 0.14.0 2020-11-22 [1] CRAN (R 4.1.0)
## glue 1.5.1 2021-11-30 [1] CRAN (R 4.1.0)
## goftest 1.2-3 2021-10-07 [1] CRAN (R 4.1.0)
## gprofiler2 * 0.2.1 2021-08-23 [1] CRAN (R 4.1.0)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.1.0)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.1.0)
## haven 2.4.3 2021-08-04 [1] CRAN (R 4.1.0)
## highr 0.9 2021-04-16 [1] CRAN (R 4.1.0)
## hms 1.1.1 2021-09-26 [1] CRAN (R 4.1.0)
## htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.1.0)
## htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.1.0)
## httpuv 1.6.3 2021-09-09 [1] CRAN (R 4.1.0)
## httr 1.4.2 2020-07-20 [1] CRAN (R 4.1.0)
## ica 1.0-2 2018-05-24 [1] CRAN (R 4.1.0)
## igraph 1.2.9 2021-11-23 [1] CRAN (R 4.1.0)
## IRanges * 2.28.0 2021-10-26 [1] Bioconductor
## irlba 2.3.5 2021-12-06 [1] CRAN (R 4.1.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
## jsonlite 1.7.2 2020-12-09 [1] CRAN (R 4.1.0)
## KEGGREST 1.34.0 2021-10-26 [1] Bioconductor
## KernSmooth 2.23-20 2021-05-03 [1] CRAN (R 4.1.0)
## knitr 1.36 2021-09-29 [1] CRAN (R 4.1.0)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.1.0)
## lattice 0.20-45 2021-09-22 [1] CRAN (R 4.1.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.1.0)
## leiden 0.3.9 2021-07-27 [1] CRAN (R 4.1.0)
## lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.1.0)
## limma 3.50.0 2021-10-26 [1] Bioconductor
## listenv 0.8.0 2019-12-05 [1] CRAN (R 4.1.0)
## lmtest 0.9-39 2021-11-07 [1] CRAN (R 4.1.0)
## lubridate 1.8.0 2021-10-07 [1] CRAN (R 4.1.0)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.1.0)
## MASS 7.3-54 2021-05-03 [1] CRAN (R 4.1.0)
## Matrix 1.4-0 2021-12-08 [1] CRAN (R 4.1.0)
## MatrixGenerics * 1.6.0 2021-10-26 [1] Bioconductor
## matrixStats * 0.61.0 2021-09-17 [1] CRAN (R 4.1.0)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.1.0)
## mgcv 1.8-38 2021-10-06 [1] CRAN (R 4.1.0)
## mime 0.12 2021-09-28 [1] CRAN (R 4.1.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.0)
## modelr 0.1.8 2020-05-19 [1] CRAN (R 4.1.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
## nlme 3.1-153 2021-09-07 [1] CRAN (R 4.1.0)
## parallelly 1.29.0 2021-11-21 [1] CRAN (R 4.1.0)
## patchwork * 1.1.1 2020-12-17 [1] CRAN (R 4.1.0)
## pbapply 1.5-0 2021-09-16 [1] CRAN (R 4.1.0)
## pillar 1.6.4 2021-10-18 [1] CRAN (R 4.1.0)
## pkgbuild 1.3.0 2021-12-09 [1] CRAN (R 4.1.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
## pkgload 1.2.4 2021-11-30 [1] CRAN (R 4.1.0)
## plotly 4.10.0 2021-10-09 [1] CRAN (R 4.1.0)
## plyr 1.8.6 2020-03-03 [1] CRAN (R 4.1.0)
## png 0.1-7 2013-12-03 [1] CRAN (R 4.1.0)
## polyclip 1.10-0 2019-03-14 [1] CRAN (R 4.1.0)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.1.0)
## princurve * 2.1.6 2021-01-18 [1] CRAN (R 4.1.0)
## processx 3.5.2 2021-04-30 [1] CRAN (R 4.1.0)
## progress 1.2.2 2019-05-16 [1] CRAN (R 4.1.0)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.1.0)
## ps 1.6.0 2021-02-28 [1] CRAN (R 4.1.0)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.1.0)
## R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.1.0)
## R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.1.0)
## R.utils 2.11.0 2021-09-26 [1] CRAN (R 4.1.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
## RANN 2.6.1 2019-01-08 [1] CRAN (R 4.1.0)
## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.1.0)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 4.1.0)
## Rcpp 1.0.7 2021-07-07 [1] CRAN (R 4.1.0)
## RcppAnnoy 0.0.19 2021-07-30 [1] CRAN (R 4.1.0)
## RCurl 1.98-1.5 2021-09-17 [1] CRAN (R 4.1.0)
## readr * 2.1.1 2021-11-30 [1] CRAN (R 4.1.0)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 4.1.0)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.1.0)
## reprex 2.0.1 2021-08-05 [1] CRAN (R 4.1.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.1.0)
## reticulate 1.22 2021-09-17 [1] CRAN (R 4.1.0)
## rhdf5 * 2.38.0 2021-10-26 [1] Bioconductor
## rhdf5filters 1.6.0 2021-10-26 [1] Bioconductor
## Rhdf5lib 1.16.0 2021-10-26 [1] Bioconductor
## rlang 0.4.12 2021-10-18 [1] CRAN (R 4.1.0)
## rmarkdown 2.11 2021-09-14 [1] CRAN (R 4.1.0)
## ROCR 1.0-11 2020-05-02 [1] CRAN (R 4.1.0)
## rpart 4.1-15 2019-04-12 [1] CRAN (R 4.1.0)
## rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.1.0)
## RSpectra 0.16-0 2019-12-01 [1] CRAN (R 4.1.0)
## RSQLite 2.2.9 2021-12-06 [1] CRAN (R 4.1.0)
## rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.1.0)
## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.1.0)
## Rtsne 0.15 2018-11-10 [1] CRAN (R 4.1.0)
## rvest 1.0.2 2021-10-16 [1] CRAN (R 4.1.0)
## S4Vectors * 0.32.3 2021-11-21 [1] Bioconductor
## sass 0.4.0 2021-05-12 [1] CRAN (R 4.1.0)
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.1.0)
## scattermore 0.7 2020-11-24 [1] CRAN (R 4.1.0)
## sctransform 0.3.2 2020-12-16 [1] CRAN (R 4.1.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.0)
## Seurat * 4.0.5 2021-10-17 [1] CRAN (R 4.1.0)
## SeuratObject * 4.0.4 2021-11-23 [1] CRAN (R 4.1.0)
## SeuratWrappers * 0.3.0 2021-06-22 [1] Github (satijalab/seurat-wrappers@fba0662)
## shiny 1.7.1 2021-10-02 [1] CRAN (R 4.1.0)
## SingleCellExperiment * 1.16.0 2021-10-26 [1] Bioconductor
## slingshot * 2.2.0 2021-10-26 [1] Bioconductor
## spatstat.core 2.3-2 2021-11-26 [1] CRAN (R 4.1.0)
## spatstat.data 2.1-0 2021-03-21 [1] CRAN (R 4.1.0)
## spatstat.geom 2.3-1 2021-12-10 [1] CRAN (R 4.1.0)
## spatstat.sparse 2.0-0 2021-03-16 [1] CRAN (R 4.1.0)
## spatstat.utils 2.3-0 2021-12-12 [1] CRAN (R 4.1.0)
## stringi 1.7.6 2021-11-29 [1] CRAN (R 4.1.0)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.1.0)
## SummarizedExperiment * 1.24.0 2021-10-26 [1] Bioconductor
## survival 3.2-13 2021-08-24 [1] CRAN (R 4.1.0)
## tensor 1.5 2012-05-05 [1] CRAN (R 4.1.0)
## testthat 3.1.1 2021-12-03 [1] CRAN (R 4.1.0)
## tibble * 3.1.6 2021-11-07 [1] CRAN (R 4.1.0)
## tidyr * 1.1.4 2021-09-27 [1] CRAN (R 4.1.0)
## tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.1.0)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.1.0)
## TrajectoryUtils * 1.2.0 2021-10-26 [1] Bioconductor
## tzdb 0.2.0 2021-10-27 [1] CRAN (R 4.1.0)
## usethis 2.1.5 2021-12-09 [1] CRAN (R 4.1.0)
## utf8 1.2.2 2021-07-24 [1] CRAN (R 4.1.0)
## uwot 0.1.11 2021-12-02 [1] CRAN (R 4.1.0)
## vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.1.0)
## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.1.0)
## viridis 0.6.2 2021-10-13 [1] CRAN (R 4.1.0)
## viridisLite 0.4.0 2021-04-13 [1] CRAN (R 4.1.0)
## withr 2.4.3 2021-11-30 [1] CRAN (R 4.1.0)
## xfun 0.28 2021-11-04 [1] CRAN (R 4.1.0)
## XML 3.99-0.8 2021-09-17 [1] CRAN (R 4.1.0)
## xml2 1.3.3 2021-11-30 [1] CRAN (R 4.1.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
## XVector 0.34.0 2021-10-26 [1] Bioconductor
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.1.0)
## zlibbioc 1.40.0 2021-10-26 [1] Bioconductor
## zoo 1.8-9 2021-03-09 [1] CRAN (R 4.1.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────